# Matrix eQTL by Andrey A. Shabalin
# http://www.bios.unc.edu/research/genomic_software/Matrix_eQTL/
# 
# Be sure to use an up to date version of R (Revolution R and Matlab are faster on PC).
#
# Set working directory:
# setwd('');
#

source('Matrix_eQTL_R/Matrix_eQTL_engine.r');
# library(MatrixEQTL) 

## Settings

# Linear model to use, modelANOVA or modelLINEAR
useModel = modelLINEAR; # modelANOVA or modelLINEAR

# Genotype file name
SNP_file_name = 'Sample_Data/SNP.txt';
snps_location_file_name = 'Sample_Data/snpsloc.txt';

# Gene expression file name
expression_file_name = 'Sample_Data/GE.txt';
gene_location_file_name = 'Sample_Data/geneloc.txt';

# Covariates file name
# Set to character() for no covariates
covariates_file_name = 'Sample_Data/Covariates.txt';

# Output file name
output_file_name = 'Sample_Data/eQTL_results_R_cis.txt';

# Only associations significant at this level will be saved
pvOutputThreshold_cis = 1e-1;
pvOutputThreshold_tra = 1e-3;

# Error covariance matrix
# Set to character() for identity.
errorCovariance = character();
# errorCovariance = read.table('Sample_Data/errorCovariance.txt');

cisDist = 1e6;

## Load genotype data

snps = SlicedData$new();
snps$fileDelimiter = '\t'; # the TAB character
snps$fileOmitCharacters = 'NA'; # denote missing values;
snps$fileSkipRows = 1; # one row of column labels
snps$fileSkipColumns = 1; # one column of row labels
snps$fileSliceSize = 2000; # read file in pieces of 2,000 rows
snps$LoadFile(SNP_file_name);

## Load gene expression data

gene = SlicedData$new();
gene$fileDelimiter = '\t'; # the TAB character
gene$fileOmitCharacters = 'NA'; # denote missing values;
gene$fileSkipRows = 1; # one row of column labels
gene$fileSkipColumns = 1; # one column of row labels
gene$fileSliceSize = 2000; # read file in pieces of 2,000 rows
gene$LoadFile(expression_file_name);

## Load covariates

cvrt = SlicedData$new();
cvrt$fileDelimiter = '\t'; # the TAB character
cvrt$fileOmitCharacters = 'NA'; # denote missing values;
cvrt$fileSkipRows = 1; # one row of column labels
cvrt$fileSkipColumns = 1; # one column of row labels
cvrt$fileSliceSize = 2000; # read file in one piece
if(length(covariates_file_name)>0) {
	cvrt$LoadFile(covariates_file_name);
}

## Run the analysis
snpspos = read.table(snps_location_file_name, header = TRUE, stringsAsFactors = FALSE);
genepos = read.table(gene_location_file_name, header = TRUE, stringsAsFactors = FALSE);

Matrix_eQTL_engine_cis(snps, 
					gene, 
					cvrt, 
					output_file_name, 
					pvOutputThreshold_cis,
					pvOutputThreshold_tra,
					useModel, 
					errorCovariance, 
					verbose = TRUE, 
					snpspos, 
					genepos,
					cisDist)